/****************************************************************************
 * Copyright (C) 2009-2010 SciTouch LLC
 * 
 * This file is part of Indigo toolkit.
 * 
 * This file may be distributed and/or modified under the terms of the
 * GNU General Public License version 3 as published by the Free Software
 * Foundation and appearing in the file LICENSE.GPL included in the
 * packaging of this file.
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 ***************************************************************************/

#include "molecule/smiles_saver.h"
#include "molecule/molecule_rgroups.h"
#include "base_cpp/array.h"
#include "base_cpp/tlscont.h"
#include "base_cpp/output.h"
#include "molecule/molecule.h"
#include "molecule/molecule_stereocenters.h"
#include "molecule/molecule_decomposer.h"
#include "graph/dfs_walk.h"

SmilesSaver::SmilesSaver (Output &output) : _output(output),
TL_CP_GET(_hcount),
TL_CP_GET(_dbonds),
TL_CP_GET(_written_atoms),
TL_CP_GET(_written_bonds)
{
   vertex_ranks = 0;
   atom_atom_mapping = 0;
   ignore_hydrogens = false;
   canonize_chiralities = false;
   write_extra_info = true;
   highlighting = 0;
}

void SmilesSaver::saveMolecule (const Molecule &mol)
{
   QS_DEF(Array<_Atom>, atoms);
   QS_DEF(Array<int>, ignored_vertices);
   int i, j, k;

   _touched_cistransbonds = 0;
   _markCisTrans(mol);

   atoms.clear_resize(mol.vertexEnd());
   atoms.zerofill();
   ignored_vertices.clear_resize(mol.vertexEnd());
   ignored_vertices.zerofill();
   _written_atoms.clear();
   _written_bonds.clear();
   _written_components = 0;
   
   if (ignore_hydrogens)
      for (i = mol.vertexBegin(); i < mol.vertexEnd(); i = mol.vertexNext(i))
         if (mol.getAtom(i).pureHydrogen() && mol.getVertex(i).degree() == 1)
            if (mol.getAtom(mol.getVertex(i).neiVertex(mol.getVertex(i).neiBegin())).label != ELEM_H)
               ignored_vertices[i] = 1;

   _hcount.clear_resize(mol.vertexEnd());

   for (i = mol.vertexBegin(); i < mol.vertexEnd(); i = mol.vertexNext(i))
   {
      _hcount[i] = mol.getAtom(i).implicit_h;

      const Vertex &vertex = mol.getVertex(i);

      if (ignore_hydrogens)
         for (j = vertex.neiBegin(); j != vertex.neiEnd(); j = vertex.neiNext(j))
            if (mol.getAtom(vertex.neiVertex(j)).pureHydrogen())
               _hcount[i]++;
   }

   
   DfsWalk walk(mol);
   
   walk.ignored_vertices = ignored_vertices.ptr();
   walk.vertex_ranks = vertex_ranks;
   walk.walk();
   
   // Detect atoms that have aromatic bonds
   for (i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
   {
      const Edge &edge = mol.getEdge(i);

      if (mol.getBond(i).type == BOND_AROMATIC)
      {
         atoms[edge.beg].aromatic = true;
         atoms[edge.end].aromatic = true;
      }
   }

   if (mol.haveQueryAtoms())
      for (i = mol.vertexBegin(); i < mol.vertexEnd(); i = mol.vertexNext(i))
      {
         if (mol.getQueryAtom(i).type != 0)
            atoms[i].aromatic = false;
      }

   // detect chiral configurations
   MoleculeStereocenters &stereocenters = mol.getStereocenters();
   const Array<DfsWalk::SeqElem> &v_seq = walk.getSequence();

   for (i = stereocenters.begin(); i != stereocenters.end(); i = stereocenters.next(i))
   {
      int atom_idx, type, group, pyramid[4];

      stereocenters.get(i, atom_idx, type, group, pyramid);

      if (type < MoleculeStereocenters::ATOM_AND)
         continue;

      int implicit_h_idx = -1;

      if (pyramid[3] == -1)
         implicit_h_idx = 3;
      else for (j = 0; j < 4; j++)
         if (ignored_vertices[pyramid[j]])
         {
            implicit_h_idx = j;
            break;
         }

      int pyramid_mapping[4];
      int counter = 0;
      QS_DEF(Array<int>, used);

      used.clear_resize(mol.vertexEnd());
      used.zerofill();

      for (j = 0; j < v_seq.size(); j++)
      {
         int v_idx = v_seq[j].idx;
         int e_idx = v_seq[j].parent_edge;

         if (e_idx >= 0 && walk.edgeClosingCycle(e_idx))
            continue;

         if (v_idx == atom_idx)
         {
            if (implicit_h_idx != -1)
               pyramid_mapping[counter++] = implicit_h_idx;

            for (int l = j + 1; l < v_seq.size(); l++)
            {
               if (v_seq[l].idx == atom_idx) // closure?
               {
                  int idx = mol.getEdge(v_seq[l].parent_edge).findOtherEnd(atom_idx);

                  used[idx] = 1;

                  for (k = 0; k < 4; k++)
                     if (pyramid[k] == idx)
                     {
                        if (counter >= 4)
                           throw Error("internal: pyramid overflow");

                        pyramid_mapping[counter++] = k;
                        break;
                     }
               }
            }

         }
         else if (!used[v_idx])
         {
            used[v_idx] = 1;
            for (k = 0; k < 4; k++)
               if (pyramid[k] == v_idx)
               {
                  if (counter >= 4)
                     throw Error("internal: pyramid overflow");
                  pyramid_mapping[counter++] = k;
                  break;
               }
         }
      }

      if (counter == 4)
      {
         // move the 'from' atom to the end
         counter = pyramid_mapping[0];
         pyramid_mapping[0] = pyramid_mapping[1];
         pyramid_mapping[1] = pyramid_mapping[2];
         pyramid_mapping[2] = pyramid_mapping[3];
         pyramid_mapping[3] = counter;
      }
      else
      {
         if (counter < 3 || _hcount[atom_idx] == 1)
            throw Error("cannot calculate chirality");
         pyramid_mapping[3] = 3;
      }

      if (MoleculeStereocenters::isPyramidMappingRigid(pyramid_mapping))
         atoms[atom_idx].chirality = 1;
      else
         atoms[atom_idx].chirality = 2;
   }

   for (i = 0; i < v_seq.size(); i++)
   {
      int v_idx = v_seq[i].idx;
      int e_idx = v_seq[i].parent_edge;

      if (e_idx == -1 || !walk.edgeClosingCycle(e_idx))
         _written_atoms.push(v_idx);

      if (e_idx != -1)
         _written_bonds.push(e_idx);
   }

   if (canonize_chiralities)
   {
      int i, j;
      QS_DEF(Array<int>, marked);
      QS_DEF(Array<int>, ids);
      const MoleculeStereocenters &stereocenters = mol.getStereocenters();

      marked.clear_resize(mol.vertexEnd());
      marked.zerofill();

      for (i = 0; i < _written_atoms.size(); i++)
      {
         if (marked[i])
            continue;

         int idx = _written_atoms[i];
         
         if (atoms[idx].chirality == 0)
            continue;

         int type = stereocenters.getType(idx);

         if (type != MoleculeStereocenters::ATOM_AND && type != MoleculeStereocenters::ATOM_OR)
            continue;

         ids.clear();
         ids.push(idx);

         int group = stereocenters.getGroup(idx);

         for (j = i + 1; j < _written_atoms.size(); j++)
         {
            if (marked[j])
               continue;

            int idx2 = _written_atoms[j];

            if (atoms[idx2].chirality == 0)
               continue;

            int type2 = stereocenters.getType(idx2);
            int group2 = stereocenters.getGroup(idx2);

            if (type2 == type && group2 == group)
            {
               ids.push(idx2);
               marked[j] = 1;
            }
         }

         if (atoms[ids[0]].chirality == 1)
            for (j = 0; j < ids.size(); j++)
               atoms[ids[j]].chirality = 3 - atoms[ids[j]].chirality;
      }
   }

   // write the SMILES itself
   
   // cycle_numbers[i] == -1 means that the number is available
   // cycle_numbers[i] == n means that the number is used by vertex n
   QS_DEF(Array<int>, cycle_numbers);

   cycle_numbers.clear();
   cycle_numbers.push(0); // never used

   bool first_component = true;
   
   for (i = 0; i < v_seq.size(); i++)
   {
      int v_idx = v_seq[i].idx;
      int e_idx = v_seq[i].parent_edge;
      int v_prev_idx = v_seq[i].parent_vertex;
      bool write_atom = true;

      if (v_prev_idx >= 0)
      {
         if (walk.numBranches(v_prev_idx) > 1)
            if (atoms[v_prev_idx].branch_cnt > 0 && atoms[v_prev_idx].paren_written)
               _output.writeChar(')');

         int opening_cycles = walk.numOpeningCycles(e_idx);
         
         for (j = 0; j < opening_cycles; j++)
         {
            for (k = 1; k < cycle_numbers.size(); k++)
               if (cycle_numbers[k] == -1)
                  break;
            if (k == cycle_numbers.size())
               cycle_numbers.push(v_prev_idx);
            else
               cycle_numbers[k] = v_prev_idx;

            _writeCycleNumber(k);
         }

         if (v_prev_idx >= 0)
         {
            int branches = walk.numBranches(v_prev_idx);
            
            if (branches > 1)
               if (atoms[v_prev_idx].branch_cnt < branches - 1)
               {
                  if (walk.edgeClosingCycle(e_idx))
                     atoms[v_prev_idx].paren_written = false;
                  else
                  {
                     _output.writeChar('(');
                     atoms[v_prev_idx].paren_written = true;
                  }
               }

            atoms[v_prev_idx].branch_cnt++;

            if (atoms[v_prev_idx].branch_cnt > branches)
               throw Error("unexpected branch");
         }

         const Bond &bond = mol.getBond(e_idx);
         const Edge &edge = mol.getEdge(e_idx);
         bool bond_written = true;

         int dir = 0;

         if (bond.type == BOND_SINGLE)
            dir = _calcBondDirection(mol, e_idx, v_prev_idx);

         if ((dir == 1 && v_idx == edge.end) || (dir == 2 && v_idx == edge.beg))
            _output.writeChar('/');
         else if ((dir == 2 && v_idx == edge.end) || (dir == 1 && v_idx == edge.beg))
            _output.writeChar('\\');
         else if (bond.type == BOND_DOUBLE)
            _output.writeChar('=');
         else if (bond.type == BOND_TRIPLE)
            _output.writeChar('#');
         else if (bond.type == BOND_AROMATIC &&
                 (!atoms[edge.beg].aromatic || !atoms[edge.end].aromatic || !mol.bondInRing(e_idx)))
            _output.writeChar(':');
         else if (bond.type == BOND_SINGLE && atoms[edge.beg].aromatic && atoms[edge.end].aromatic)
            _output.writeChar('-');
         else
            bond_written = false;

         if (walk.edgeClosingCycle(e_idx))
         {
            for (j = 1; j < cycle_numbers.size(); j++)
               if (cycle_numbers[j] == v_idx)
                  break;

            if (j == cycle_numbers.size())
               throw Error("cycle number not found");

            _writeCycleNumber(j);

            cycle_numbers[j] = -1;
            write_atom = false;
         }
      }
      else
      {
         if (!first_component)
            _output.writeChar('.');
         first_component = false;
         _written_components++;
      }
      if (write_atom)
         _writeAtom(mol, v_idx, atoms[v_idx].aromatic, atoms[v_idx].chirality);
   }

   if (write_extra_info)
   {
      _comma = false;
      _writeStereogroups(mol, atoms);
      _writeRadicals(mol);
      _writePseudoAtoms(mol);
      _writeHighlighting();

      if (_comma)
         _output.writeChar('|');
   }
}

void SmilesSaver::_writeCycleNumber (int n) const
{
   if (n > 0 && n < 10)
      _output.printf("%d", n);
   else if (n >= 10 && n < 100)
      _output.printf("%%%2d", n);
   else if (n >= 100 && n < 1000)
      _output.printf("%%%%%3d", n);
   else 
      throw Error("bad cycle number: %d", n);
}

void SmilesSaver::_writeAtom (const Molecule &mol, int idx, bool aromatic, int chirality) const
{
   const Atom &atom = mol.getAtom(idx);
   const QueryAtom *query_atom = 0;
   int i;
   bool need_brackets = false;
   int hydro = -1;
   int aam = 0;

   if (mol.haveQueryAtoms())
   {
      query_atom = &mol.getQueryAtom(idx);

      if (query_atom->type == QUERY_ATOM_RGROUP)
      {
         if (mol.getRGroups()->isRGroupAtom(idx))
         {
            const Array<int> &rg = mol.getRGroups()->getSiteRGroups(idx);

            if (rg.size() != 1)
               throw Error("rgroup count %d", rg.size());
            
            _output.printf("[&%d]", rg[0] + 1);
         }
         else
            _output.printf("[&%d]", 1);

         return;
      }

      if (query_atom->type == QUERY_ATOM_A)
      {
         _output.printf("[*]");
         return;
      }
   }

   if (mol.isPseudoAtom(idx))
   {
      _output.printf("*");
      return;
   }

   if (atom_atom_mapping != 0)
      aam = atom_atom_mapping[idx];

   if (atom.label != ELEM_C && atom.label != ELEM_P &&
       atom.label != ELEM_N && atom.label != ELEM_S &&
       atom.label != ELEM_O && atom.label != ELEM_Cl &&
       atom.label != ELEM_F && atom.label != ELEM_Br &&
       atom.label != ELEM_B && atom.label != ELEM_I)
      need_brackets = true;

   if (atom.explicit_valence || atom.radical != 0 || chirality > 0 ||
       (aromatic && atom.label != ELEM_C && atom.label != ELEM_O) ||
       (aromatic && atom.label == ELEM_C && mol.getVertex(idx).degree() < 3 && _hcount[idx] == 0))
      hydro = _hcount[idx];

   if (chirality || atom.charge != 0 || atom.isotope > 0 || hydro >= 0 || aam > 0)
      need_brackets = true;

   if (need_brackets)
   {
      if (hydro == -1)
         hydro = _hcount[idx];
      _output.writeChar('[');
   }

   if (atom.isotope > 0)
      _output.printf("%d", atom.isotope);

   const char *elem = Element::toString(atom.label);

   if (aromatic)
   {
      for (i = 0; i < (int)strlen(elem); i++)
         _output.printf("%c", tolower(elem[i]));
   }
   else
      _output.printf("%s", elem);

   if (chirality > 0)
   {
      if (chirality == 1)
         _output.printf("@");
      else // chirality == 2
         _output.printf("@@");
      
      if (atom.implicit_h > 1)
         throw Error("%d implicit H near stereocenter", atom.implicit_h);
   }

   if (hydro > 1)
      _output.printf("H%d", hydro);
   else if (hydro == 1)
      _output.printf("H");

   if (atom.charge > 1)
      _output.printf("+%d", atom.charge);
   else if (atom.charge < -1)
      _output.printf("-%d", -atom.charge);
   else if (atom.charge == 1)
      _output.printf("+");
   else if (atom.charge == -1)
      _output.printf("-");

   if (aam > 0)
      _output.printf(":%d", aam);

   if (need_brackets)
      _output.writeChar(']');

   if (mol.getRGroupFragment() != 0)
   {
      for (i = 0; i < 2; i++)
      {
         int j;

         for (j = 0; mol.getRGroupFragment()->getAttachmentPoint(i, j) != -1; j++)
            if (idx == mol.getRGroupFragment()->getAttachmentPoint(i, j))
            {
               _output.printf("([*])");
               break;
            }

         if (mol.getRGroupFragment()->getAttachmentPoint(i, j) != -1)
            break;
      }
   }
}

void SmilesSaver::_markCisTrans (const Molecule &mol)
{
   int i, j;

   _dbonds.clear_resize(mol.edgeEnd());

   for (i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
   {
      _dbonds[i].ctbond_beg = -1;
      _dbonds[i].ctbond_end = -1;
      _dbonds[i].saved = 0;
   }

   if (!mol.cis_trans.exists())
      return;

   for (i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
      if (mol.cis_trans.getParity(i) != 0 && !mol.bondInRing(i))
      {
         const Edge &edge = mol.getEdge(i);
         const Vertex &beg = mol.getVertex(edge.beg);
         const Vertex &end = mol.getVertex(edge.end);
         bool arom_fail_beg = true, arom_fail_end = true;

         for (j = beg.neiBegin(); j != beg.neiEnd(); j = beg.neiNext(j))
         {
            int idx = beg.neiEdge(j);

            if (idx != i && mol.getBond(idx).type == BOND_SINGLE)
               arom_fail_beg = false;
         }

         for (j = end.neiBegin(); j != end.neiEnd(); j = end.neiNext(j))
         {
            int idx = end.neiEdge(j);

            if (idx != i && mol.getBond(idx).type == BOND_SINGLE)
               arom_fail_end = false;
         }

         if (arom_fail_beg || arom_fail_end)
            continue;


         for (j = beg.neiBegin(); j != beg.neiEnd(); j = beg.neiNext(j))
         {
            int idx = beg.neiEdge(j);

            if (idx != i)
            {
               if (mol.getEdge(idx).beg == edge.beg)
                  _dbonds[idx].ctbond_beg = i;
               else
                  _dbonds[idx].ctbond_end = i;
            }
         }

         for (j = end.neiBegin(); j != end.neiEnd(); j = end.neiNext(j))
         {
            int idx = end.neiEdge(j);

            if (idx != i)
            {
               if (mol.getEdge(idx).beg == edge.end)
                  _dbonds[idx].ctbond_beg = i;
               else
                  _dbonds[idx].ctbond_end = i;
            }
         }
      }
}

bool SmilesSaver::_updateSideBonds (const Molecule &mol, int bond_idx)
{
   const Edge &edge = mol.getEdge(bond_idx);
   int subst[4];
   mol.cis_trans.getSubstituents_All(mol, bond_idx, subst);
   int parity = mol.cis_trans.getParity(bond_idx);

   int sidebonds[4] = {-1, -1, -1, -1};

   sidebonds[0] = mol.findEdgeIndex(subst[0], edge.beg);
   if (subst[1] != -1)
      sidebonds[1] = mol.findEdgeIndex(subst[1], edge.beg);

   sidebonds[2] = mol.findEdgeIndex(subst[2], edge.end);
   if (subst[3] != -1)
      sidebonds[3] = mol.findEdgeIndex(subst[3], edge.end);

   int n1 = 0, n2 = 0, n3 = 0, n4 = 0;

   if (_dbonds[sidebonds[0]].saved != 0)
   {
      if ((_dbonds[sidebonds[0]].saved == 1 && mol.getEdge(sidebonds[0]).beg == edge.beg) ||
          (_dbonds[sidebonds[0]].saved == 2 && mol.getEdge(sidebonds[0]).end == edge.beg))
         n1++;
      else
         n2++;
   }
   if (sidebonds[1] != -1 && _dbonds[sidebonds[1]].saved != 0)
   {
      if ((_dbonds[sidebonds[1]].saved == 2 && mol.getEdge(sidebonds[1]).beg == edge.beg) ||
          (_dbonds[sidebonds[1]].saved == 1 && mol.getEdge(sidebonds[1]).end == edge.beg))
         n1++;
      else
         n2++;
   }
   if (_dbonds[sidebonds[2]].saved != 0)
   {
      if ((_dbonds[sidebonds[2]].saved == 1 && mol.getEdge(sidebonds[2]).beg == edge.end) ||
          (_dbonds[sidebonds[2]].saved == 2 && mol.getEdge(sidebonds[2]).end == edge.end))
         n3++;
      else
         n4++;
   }
   if (sidebonds[3] != -1 && _dbonds[sidebonds[3]].saved != 0)
   {
      if ((_dbonds[sidebonds[3]].saved == 2 && mol.getEdge(sidebonds[3]).beg == edge.end) ||
          (_dbonds[sidebonds[3]].saved == 1 && mol.getEdge(sidebonds[3]).end == edge.end))
         n3++;
      else
         n4++;
   }

   if (parity == MoleculeCisTrans::CIS)
   {
      n1 += n3;
      n2 += n4;
   }
   else
   {
      n1 += n4;
      n2 += n3;
   }

   if (n1 > 0 && n2 > 0)
      throw Error("incompatible cis-trans configuration");

   if (n1 == 0 && n2 == 0)
      return false;

   if (n1 > 0)
   {
      _dbonds[sidebonds[0]].saved =
         (mol.getEdge(sidebonds[0]).beg == edge.beg) ? 1 : 2;
      if (sidebonds[1] != -1)
         _dbonds[sidebonds[1]].saved =
            (mol.getEdge(sidebonds[1]).beg == edge.beg) ? 2 : 1;

      _dbonds[sidebonds[2]].saved =
         ((mol.getEdge(sidebonds[2]).beg == edge.end) == (parity == MoleculeCisTrans::CIS)) ? 1 : 2;
      if (sidebonds[3] != -1)
         _dbonds[sidebonds[3]].saved =
            ((mol.getEdge(sidebonds[3]).beg == edge.end) == (parity == MoleculeCisTrans::CIS)) ? 2 : 1;
   }
   if (n2 > 0)
   {
      _dbonds[sidebonds[0]].saved =
         (mol.getEdge(sidebonds[0]).beg == edge.beg) ? 2 : 1;
      if (sidebonds[1] != -1)
         _dbonds[sidebonds[1]].saved =
            (mol.getEdge(sidebonds[1]).beg == edge.beg) ? 1 : 2;

      _dbonds[sidebonds[2]].saved =
         ((mol.getEdge(sidebonds[2]).beg == edge.end) == (parity == MoleculeCisTrans::CIS)) ? 2 : 1;
      if (sidebonds[3] != -1)
         _dbonds[sidebonds[3]].saved =
            ((mol.getEdge(sidebonds[3]).beg == edge.end) == (parity == MoleculeCisTrans::CIS)) ? 1 : 2;
   }

   return true;
}

int SmilesSaver::_calcBondDirection (const Molecule &mol, int idx, int vprev)
{
   int i, ntouched;

   if (_dbonds[idx].ctbond_beg == -1 && _dbonds[idx].ctbond_end == -1)
      return 0;

   if (mol.getBond(idx).type != BOND_SINGLE)
      throw Error("internal: directed bond type %d", mol.getBond(idx).type);

   while (1)
   {
      ntouched = 0;
      for (i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
         if (mol.cis_trans.getParity(i) != 0 && !mol.bondInRing(i))
         {
            if (_updateSideBonds(mol, i))
               ntouched++;
         }
      if (ntouched == _touched_cistransbonds)
         break;
      _touched_cistransbonds = ntouched;
   }

   if (_dbonds[idx].saved == 0)
   {
      if (vprev == mol.getEdge(idx).beg)
         _dbonds[idx].saved = 1;
      else
         _dbonds[idx].saved = 2;
   }

   return _dbonds[idx].saved;
}

void SmilesSaver::_writeStereogroups (const Molecule &mol, const Array<_Atom> &atoms)
{
   MoleculeStereocenters &stereocenters = mol.getStereocenters();
   int i, j;
   int single_and_group = -1;

   for (i = stereocenters.begin(); i != stereocenters.end(); i = stereocenters.next(i))
   {
      int idx, type, group;

      stereocenters.get(i, idx, type, group, 0);

      if (type < MoleculeStereocenters::ATOM_ANY)
         continue;
      if (type != MoleculeStereocenters::ATOM_AND)
         break;
      if (single_and_group == -1)
         single_and_group = group;
      else if (single_and_group != group)
         break;
   }

   if (i == stereocenters.end())
      return;

   int and_group_idx = 1;
   int or_group_idx = 1;
   
   QS_DEF(Array<int>, marked);

   marked.clear_resize(_written_atoms.size());
   marked.zerofill();

   for (i = 0; i < _written_atoms.size(); i++)
   {
      if (marked[i])
         continue;

      int type = stereocenters.getType(_written_atoms[i]);

      if (type > 0)
      {
         if (_comma)
            _output.writeChar(',');
         else
         {
            _output.writeString(" |");
            _comma = true;
         }
      }

      if (type == MoleculeStereocenters::ATOM_ANY)
      {
         _output.printf("w:%d", i);

         for (j = i + 1; j < _written_atoms.size(); j++)
            if (stereocenters.getType(_written_atoms[j]) == MoleculeStereocenters::ATOM_ANY)
            {
               marked[j] = 1;
               _output.printf(",%d", j);
            }
      }
      else if (type == MoleculeStereocenters::ATOM_ABS)
      {
         _output.printf("a:%d", i);

         for (j = i + 1; j < _written_atoms.size(); j++)
            if (stereocenters.getType(_written_atoms[j]) == MoleculeStereocenters::ATOM_ABS)
            {
               marked[j] = 1;
               _output.printf(",%d", j);
            }
      }
      else if (type == MoleculeStereocenters::ATOM_AND)
      {
         int group = stereocenters.getGroup(_written_atoms[i]);
         
         _output.printf("&%d:%d", and_group_idx++, i);
         for (j = i + 1; j < _written_atoms.size(); j++)
            if (stereocenters.getType(_written_atoms[j]) == MoleculeStereocenters::ATOM_AND &&
                stereocenters.getGroup(_written_atoms[j]) == group)
            {
               marked[j] = 1;
               _output.printf(",%d", j);
            }
      }
      else if (type == MoleculeStereocenters::ATOM_OR)
      {
         int group = stereocenters.getGroup(_written_atoms[i]);

         _output.printf("o%d:%d", or_group_idx++, i);
         for (j = i + 1; j < _written_atoms.size(); j++)
            if (stereocenters.getType(_written_atoms[j]) == MoleculeStereocenters::ATOM_OR &&
                stereocenters.getGroup(_written_atoms[j]) == group)
            {
               marked[j] = 1;
               _output.printf(",%d", j);
            }
      }
   }
}

void SmilesSaver::_writeRadicals (const Molecule &mol)
{
   QS_DEF(Array<int>, marked);
   int i, j;

   marked.clear_resize(_written_atoms.size());
   marked.zerofill();

   for (i = 0; i < _written_atoms.size(); i++)
   {
      if (marked[i])
         continue;

      int radical = mol.getAtom(_written_atoms[i]).radical;

      if (radical == 0)
         continue;

      if (_comma)
         _output.writeChar(',');
      else
      {
         _output.writeString(" |");
         _comma = true;
      }
    
      if (radical == RADICAL_SINGLET)
         _output.writeString("^3:");
      else if (radical == RADICAL_DOUPLET)
         _output.writeString("^1:");
      else // RADICAL_TRIPLET
         _output.writeString("^4:");

      _output.printf("%d", i);
      
      for (j = i + 1; j < _written_atoms.size(); j++)
         if (mol.getAtom(_written_atoms[j]).radical == radical)
         {
            marked[j] = 1;
            _output.printf(",%d", j);
         }
   }
}

void SmilesSaver::_writePseudoAtoms (const Molecule &mol)
{
   int i;
   
   for (i = 0; i < _written_atoms.size(); i++)
   {
      if (mol.isPseudoAtom(_written_atoms[i]))
         break;
   }

   if (i == _written_atoms.size())
      return;

   if (_comma)
      _output.writeChar(',');
   else
   {
      _output.writeString(" |");
      _comma = true;
   }

   _output.writeChar('$');

   for (i = 0; i < _written_atoms.size(); i++)
   {
      if (i > 0)
         _output.writeChar(';');

      if (mol.isPseudoAtom(_written_atoms[i]))
         writePseudoAtom(mol.getPseudoAtom(_written_atoms[i]), _output);
   }

   _output.writeChar('$');
}

void SmilesSaver::writePseudoAtom (const char *label, Output &out)
{
   if (*label == 0)
      throw Error("empty pseudo-atom");

   do
   {
      if (isspace(*label))
         throw Error("character 0x%x is not allowed inside pseudo-atom", *label);
      if (*label == '$' || *label == ';')
         throw Error("'%c' not allowed inside pseudo-atom", *label);

      out.writeChar(*label);
   } while (*(++label) != 0);
}

void SmilesSaver::_writeHighlighting ()
{
   if (highlighting == 0)
      return;

   int i;

   bool ha = false;

   for (i = 0; i < _written_atoms.size(); i++)
   {
      if (highlighting->hasVertex(_written_atoms[i]))
      {
         if (ha)
            _output.writeChar(',');
         else
         {
            if (_comma)
               _output.writeChar(',');
            else
            {
               _output.writeString(" |");
               _comma = true;
            }
            _output.writeString("ha:");
            ha = true;
         }

         _output.printf("%d", i);
      }
   }

   bool hb = false;

   for (i = 0; i < _written_bonds.size(); i++)
   {
      if (highlighting->hasEdge(_written_bonds[i]))
      {
         if (hb)
            _output.writeChar(',');
         else
         {
            if (_comma)
               _output.writeChar(',');
            else
            {
               _output.writeString(" |");
               _comma = true;
            }
            _output.writeString("hb:");
            hb = true;
         }

         _output.printf("%d", i);
      }
   }
}


int SmilesSaver::writtenComponents ()
{
   return _written_components;
}

const Array<int> & SmilesSaver::writtenAtoms ()
{
   return _written_atoms;
}

const Array<int> & SmilesSaver::writtenBonds ()
{
   return _written_bonds;
}
